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ABSTRACT 

The composition of Ultra High Energy Gosmic Rays (UHEGR) is measured 
with the High Resolution Fly's Eye cosmic ray observatory (HiRes) data using 
the ^max technique. Data were collected in stereo between 1999 November and 
2001 September. The data are reconstructed with well-determined geometry. 
Measurements of the atmospheric transmission are incorporated in the recon- 
struction. The detector resolution is found to be 30 g cm~^ in Xmax and 13% in 
Energy. The X^ax elongation rate between 10^^'° eV and 10^^'^ eV is measured to 
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be 54.5 ± 6.5 (stat) ± 4.5 (sys) g cm~^ per decade. This is compared to predic- 
tions using the QGSJetOl and SIBYLL 2.1 hadronic interaction models for both 
protons and iron nuclei. CORSIKA-generated Extensive Air Showers (EAS) are 
incorporated directly into a detailed detector Monte Carlo program. The elon- 
gation rate and the Xmax distribution widths are consistent with a constant or 
slowly changing and predominantly light composition. A simple model contain- 
ing only protons and iron nuclei is compared to QGSJet and SIBYLL. The best 
agreement between the model and the data is at 80% protons for QGSJet and 
60% protons for SIBYLL. 

Subject headings: cosmic rays — acceleration of particles — large-scale structure 
of universe 



1. INTRODUCTION 

The cosmic ray (CR) spectrum follows a power law which exhibits several interesting 
features. A break (known as the knee) between a power law with index -2.7 to index -3.0 
occurs near 10^^ eV. A similar break (the second knee) from index -3.0 to -3.3 has been 
reported to occur near lO^^-^ eV (Abu-Zayyad et al. 2001a; Bird et al. 1993, 1995; Pravdin 
et al. 1999, 2003; Nagano et al. 1984, 1991), followed by a rise to an index of -2.7 (the ankle) 
near 10^^'® eV. The spectrum appears to continue until 10^^'® eV. 

Changes in spectral index at various energies could be produced by a gradual decrease 
in efficiency with energy of a source (e.g., galactic supernovae), leakage out of an area 
of magnetic confinement (the galactic leaky box model), the appearance of flux from new 
sources that begin to dominate at higher energies (extra-galactic cosmic rays), and thresholds 
of inelastic interactions between CR protons and the cosmic microwave background radiation 
(e+e~ and pion production, (the GZK effect (Greisen 1966; Zatespin & Kuz'min 1966))). 
It has also been suggested that the hadronic interaction of the CR with the atmosphere 
undergoes changes above certain energies which would affect the measured energy of the CR 
and hence produce an apparent change in the power law index (Kazanas & Nicolaidis 2003). 
A knowledge of CR composition as a function of energy would be invaluable in sorting out 
these effects. 

High energy cosmic rays have long been known to be charged nuclei (Clay 1927), but 
determining the chemical composition at energies greater than 10^^ cV is especially difficult 
because of the low flux. The High Resolution Fly's Eye cosmic ray observatory has a large 
enough aperture (3x10^ to 5x10^ km^ str) but it observes the extensive air shower (EAS) 
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produced by the particle rather than detecting the primary itself, and thus must use an 
indirect method to study the composition. 

The distribution of positions of shower maxima Q^max) in the atmosphere has been 
shown to be sensitive to the composition of cosmic rays (Heitler 1944). It is well known 
that for any particular species of nucleus, the position of shower maximum will deepen with 
increasing energy as the logarithm of the energy. The slope, d{Xmax) / d{logE) , is known as 
the elongation rate (Linsley 1977). 

While the details are dependent on the hadronic model assumed, all modern hadronic 
models give approximately the same elongation rate (between 50 and 60 g cm~^ per decade 
of energy, independent of particle species) and agree within about 25 g cm~^ on the absolute 
position of the average shower yi,nax at a given energy for a given species. The sensitivity of 
the Xfnax method to composition comes from the fact that the mean X^ax for iron and protons 
is different by about 75 g cm~^, independent of hadronic model with protons producing 
deeper showers with larger fluctuations. A change in the composition from heavy to light 
would then result in a larger elongation rate than 50-60 g cm~^ per decade, and a change 
from light to heavy would lead to a lower and even negative elongation rate. 

Previous experiments (stereo Fly's Eye (Bird et al. 1993), HiRes prototype-MIA (Abu- 
Zayyad et al. 2000)) have shown evidence for an elongation rate of 80-90 g cm~^ g cm~^ 
per decade in the energy rage from 10^^ to 10^^"^ eV. No significant information from air- 
fluorescence experiments has been hitherto available on the behavior of the elongation rate 
near 10^^ eV and above. 

The general dependence of 'Kmax on energy can be seen in a simple branching model 
where Nj^ax ^ and X^ax ^ ln(£'o) (Gaisser 1990; Heitler 1944). In this model, if the 
primary particle is a nucleus the shower is assumed to be a superposition of subshowers each 
initiated by one of the A independent nucleons. The primary energy must be divided among 
the A constituents, so in this case X^ax oc \n.{A/Eo). A more complete discussion leads to 
Linsley's expression for the elongation rate, ct, (Gaisser et al. 1979; Linsley 1977a), which is 



a = (1 - B)KX 



1 



d{log{A)) 



(1) 



d{logEo] 

where K is a constant, A is the collision length, and B expresses the dependence of a on the 
hadron-air nucleus interactions. It includes both the energy dependence of the cross-section 
(and thus A) and the energy dependence of the multiplicity and inelasticity (Heck 2001; 
Linsley & Watson 1981; Feynman 1969). 



The technique for extracting the cosmic ray composition used in this paper reduces to 
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comparing the ^max distribution of the data, after appropriate cuts that guarantee good 
resolution in this variable, to simulated data, generated with either a proton or iron parent 
particle. The simulated data is the result of a detailed detector Monte Carlo and includes 
all the reconstruction uncertainties. 



2. THE HiRes DETECTOR 

The High Resolution Fly's Eye Cosmic Ray Observatory (HiRes) has two sites on the 
U.S. Army's Dugway Proving Ground in the West Desert of Utah, about 90 miles from Salt 
Lake City. The first site, HiRes- 1, is on Little Granite Mountain, the site of the original 
Fly's Eye detector. HiRes-2 is 12.6 km to the southwest on Camel's Back Ridge. Dugway 
was chosen for its clean atmosphere and low light pollution. Each site is on a hill, above 
the bulk of most haze in the atmosphere. The two sites gather data independently, and the 
data can be analyzed from each site in monocular mode or together in stereo mode. Results 
from a monocular analysis of the spectrum have been published (Abbasi et al. 2004, a). 
The accurate determination of the shower geometry is an essential first step to determining 
the shower profile, which gives the energy of the primary UHECR and the ^max of the 
shower. This work will present results from stereo analysis, which has an obvious advantage 
in determining the shower geometry and hence better resolution in energy and X^ax- 

HiRes is the realization of an extension of the method pioneered by the original Fly's 
Eye experiment (Baltrusaitis et al. 1985). As a UHECR-initiated Extensive Air Shower 
(EAS) propagates through the atmosphere, the charged particles excite nitrogen molecules, 
which fluoresce. The fluorescence yield and its spectrum have both been measured (Bunner 
1967; Kakimoto et al. 1996). The fluorescence yield is about five photons per particle per 
meter. The photons, mostly of wavelength 300-400 nm, propagate isotropically from the 
shower core, with the number of photons coming from a slice of the shower proportional 
to the number of charged particles in that slice. In addition to air fluorescence, Cerenkov 
hght is produced by shower particles. The Cerenkov cone in air is 2.5° and the resultant 
light angular distribution, reflecting the multiple scattering of electrons, is strongly beamed 
forward along the EAS axis. Some of this light, however, will scatter at larger angles into 
the detector and needs to be subtracted from the total signal. 

Air fluorescence detectors gather the photons and focus them onto arrays of photo- 
multipher tubes (PMTs) . Each PMT views a solid angle of the sky, and a composite view 
is recovered when the images from the constituent mirrors are combined. 

The HiRes Prototype detector, results from which are also summarized in this paper. 
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was located on Little Granite Mountain (Abu-Zayyad et al. 2000a). The prototype detector 
viewed from 3° to 70° in elevation and had an azimuthal coverage that overlooks the Michigan 
Muon Array (MIA) 3.4 km away. Coincident hybrid data were collected from 1993 to 1996. 
Results on composition near 10^^ eV were published in (Abu-Zayyad et al. 2000). 

In 1997, the aperture for events at the highest energies was optimized by redeploying 
the telescopes in a ring with full 360° azimuthal coverage and 3° to 17° elevational coverage 
(Abbasi et al. 2004a). HiRes-1 began taking data in May of 1997, and the ring was complete 
in March of 1998. 

HiRes-2 has two rings, giving complete azimuthal coverage with an elevation coverage 
from 3° to 31° (Boyer et al. 2002). This analysis uses data from November of 1999, which 
is when HiRes-2 was completely operational, to September of 2001. 

2.1. Survey and Calibration 

The basic HiRes element is a telescope consisting of a mirror of area 5.2 m^ with an 
associated cluster of 256 PMTs and data-acquisition (DAQ) electronics. The PMT cluster is 
placed at a distance from the mirror which optimizes the spot size. (Simpson 2001). Taking 
into account the cluster obscuration, the effective area of each mirror is 3.72 m^ (Abbasi et 
al. 2004a). The pointing directions of the individual telescopes were surveyed at installation 
and are checked periodically by observing signals produced by stars (Bergman et al. 2001; 
Sadowski et al. 2002). 

The primary tool for the calibration of PMT sensitivity is a xenon flash lamp mounted 
in a portable housing that can be moved from telescope to telescope. It is placed in the 
center of the mirror, illuminating the cluster directly. The lamp's output has been measured 
to be stable to within 1/3% flash-to-flash and within 2% over the course of a night (Jones 
et al. 1999, 2001). 

To monitor the PMT response on a nightly basis, each site has a frequency-tripled YAG 
laser which delivers light at 355 nm to each PMT cluster via quartz optical fibers. One fiber 
goes to the center of each mirror, and one to each side of each PMT cluster. The mirror- 
mounted fibers illuminate the cluster directly allowing monitoring of tube response, while 
the cluster-mounted fibers illuminate the mirrors so mirror refiectivity can be tracked (Girard 
et al. 2000; Archbold et al. 2001). 
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3. 



ATMOSPHERIC CORRECTIONS 



Photons from an EAS travel through the atmosphere to reach the telescopes, and un- 
derstanding light transmission through the atmosphere is vital. The atmosphere can be 
considered a mix of molecules and aerosols. Rayleigh scattering from molecules is well un- 
derstood, as is the atmospheric density profile. Scattering from aerosols varies with the 
aerosol content of the air and must be measured. We characterize the aerosol scattering 
by the horizontal aerosol extinction length, L", the vertical aerosol optical depth, VAGD, 
and by the angular dependence of the scattering cross section (the phase function). Lasers 
with steerable beams, located at both HiRes-1 and HiRes-2 sites, are used to sweep the 
HiRes aperture and determine the horizontal extinction length, VAOD and phase function 
(Roberts et al. 2001; Wiencke et al. 2001). A database of the hourly parameterizations of 
the atmosphere is used in Monte Carlo and for the reconstruction of every event for which 
it is available (Reil 2002). For the final data set, about three-quarters of the events were re- 
constructed with atmospheric parameters from the hourly database and the remainder were 
reconstructed with the average values (see Section 4.6). 

The laser tracks also give an indication of clouds in the aperture. The track of a laser 
hitting a cloud will mushroom out, giving a shorter, wider signature. 

Other tools for understanding the clarity of the atmosphere include operator observa- 
tions, infrared cloud monitors, and xenon flashers (Wiencke et al. 1999). One inclined and 
ten vertical xenon flashers located between the two HiRes sites are fired every ten minutes 
and give important qualitative information about the atmosphere between the sites. 



The initial steps in the data processing chain are documented by Reil (2002) and Abu- 
Zayyad (2000b). First, calibration information is applied to the raw data. The relative 
timing information from each mirror is then converted to an absolute time as determined by 
GPS (Wilkinson 1998). The individual mirror triggers are matched to form multiple- mirror 

events, and the multiple-mirror events from each site are time-matched to build stereo events. 
To separate noise events from track- like events, a Rayleigh filter is employed (Song 2001). 
The probability that an event was created by random noise is required to be less than 0.1%. 



4. 



DATA ANALYSIS 



4.1. 



Reconstruction 
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4.2. Geometry 

The axis of the EAS and the position of the detector uniquely define the Shower-Detector 
Plane (SDP), as illustrated in Fig. 1. The location and pointing direction of each PMT 
cluster have been measured (Bergman et al. 2001; Sadowski et al. 2002), so the pointing 
direction n, of each PMT is known and the SDP is easily found. 

Once the SDP is known for each site, the intersection of the planes gives the direction 
and location of the EAS. The next step in the reconstruction is to calculate the shower 
development profile. However, light arriving at the detector is collected by discrete PMTs, 
each of which covers about 1° x 1° of the sky. The signal from a longitudinal segment of the 
EAS is thus necessarily spfit among many PMTs. For profile fitting, the signal must be re- 
combined into bins that correspond to the longitudinal segments of the EAS. The re- binned 
signal, corrected for atmospheric extinction and with Cerenkov light subtracted is fit to a 
Gaisser-Hillas functional form (Eq. 2). This form has been shown to be in good agreement 
with EAS simulations (Song 2001; Heck 1998; Kalmykov et al. 1997) and with HiRes data 
(Abu-Zayyad et al. 2001). 

y ~_ ° exp 

4.3. Angulctr Binning 

A detailed description of the angular binning technique can be found in Reil (2002) and 
Archbold (2002). With the incoming photon flux divided into the angular bins, an "inverse 
Monte Carlo" method is employed to correct for acceptance. A Monte Carlo shower with 
geometry corresponding to the event in question is generated with a Gaisser-Hillas profile 
(see Eqn. 2). Each photon, including both scintillation photons and Cerenkov photons, is 
individually traced up to the same point where the fiux $ was computed from the data, giving 
^MC r^Yie atmospheric parameters described in Section 3 are used to calculate attenuation 
and scattering. If measured values of horizontal attenuation length and scale height are 
available in the data base for the hour during which the event occurred, those measured 
values are used. Otherwise, the average atmospheric values are used (Wiencke et al. 2001a). 
^max and Xmax are then allowed to vary to minimize the with respect to the measured 
fiux 



{-^max ^) 

A 



(2) 



-8- 



Xmc — 




(3) 



where j is the number of bins. For this x^, 



, (7 is given by 




+ 0"^r Ca 



(4) 



where Npe is the number of photoelectrons in the bin, C^^yy is a geometric correction factor 
for non-normal incidence, and 



where 40 (pe///sec) is the average sky noise. 

Eq. (3) can be minimized for HiRes-1 and HiRes-2 individually, or for both sites globally. 
The ^max and ^max that minimize Eq. (3) (along with the parameters Xq and A) define 
the shower profile. The total number of charged particles is obtained from the integral of 
Eq. (2), and the energy is calculated by multiplying the total number of charged particles 
by the energy deposited per charged particle. A correction is made for unobserved energy 
averaging about 10%. 



In addition to this analysis method, a reconstruction technique which takes advantage 
of the FADC timing at HiRes-2 has also been developed (Abbasi et al. 2004a). It has the 
advantage of being less dependent on the details of phototube acceptance. We adapt it to 
generate one of the quality cuts described below. 

In this method, the SDP are found as above and the EAS geometry obtained from the 
intersection of the SDP is assumed. A first estimate of the number of charged particles at 
the EAS in each HiRes-2 time bin is calculated assuming all of the photons reaching the 
detector are from air fluorescence. Note that each time bin contains the contributions of a 
number of PMTs determined by the effective optical spot size. The first guess of the number 
of charged particles at the shower obtained from the data is compared to the number of 
charged particles from a Gaisser-Hillas profile. A scan through each 'Ys.max and its associated 
most likely ^max is performed to find the best fit. Cerenkov light is then introduced, based 
on the number of charged particles from the Gaisser-Hillas fit. The Cerenkov hght is traced 
to the detector and subtracted from the signal, and the process is repeated to find a new 



(5) 



4.4. Time Binning 
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Xmaa; ^nd ^max- The iteration continues until satisfactory agreement between the predicted 
data and real data is obtained. The energy of the event is calculated from the shower profile 
exactly as in the angular binning technique. 

4.5. Stereo Reconstruction Cuts 

To insure that both detectors were working properly, only data files in which at least 
20 xenon vertical flashers events (see Section 3) are seen by both eyes are used for this 
analysis. Wc chose events of energy greater than 10^^ eV because the stereo aperture is 
rapidly decreasing below this energy. The statistical reach (defined as at least 4 events per 
bin) of this sample of data corresponds to a maximum bin at 10^'^'^ cV. Additional loose cuts 
described in Table 1 were used to remove obviously badly reconstructed events and stereo 
mismatches without biasing the data sample. For the period between November of 1999 
and September of 2001, 728 events met all of the above criteria and were subjected to the 
atmospheric cuts. 

4.6. Atmospheric Cuts 

For this analysis, any event with a corresponding VAOD measurement of larger than 0.1 
was cut. However, because of equipment down time, the atmosphere was not measured for 
every hour for which we have data. In that case the operators' log comments and measure- 
ments of laser track length vs. width (see Section 3) were searched. Periods during which 
the operators' comments indicated bad weather and/or the length vs. width of the laser 
tracks indicated that the aperture was cloudy were discarded. Of the 553 events comprising 
the final data set, 419 had atmospheric database entries which were used for reconstruction. 
The remaining 134 events had no database entry but occurred during good weather and were 
reconstructed with parameters corresponding to average atmospheric conditions. 

4.7. Quality Cuts 

All remaining events were manually scanned using an event display. The cuts listed in 
Table 1 produce a data set that contains real events during periods of operation where the 
atmosphere was well understood and the detector was working well. We use the Monte Carlo 
simulated data to define these loose cuts to insure that sufficiently precise determination of 
event energy and Xmax is had. To establish these cuts, a Monte Carlo set of 8341 events 
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was generated with a Fly's Eye Stereo spectrum for these resolution studies. The input 
composition for the Monte Carlo was nearly equal numbers of protons and iron nuclei, 
evenly divided between the QGSJet and SIBYLL hadronic interaction models to reflect the 
extremes of possible composition and hadronic interaction models. Table 1 summarizes the 
resultant cuts. 

Proper reconstruction of the shower profile for events pointed at the detector is espe- 
cially dependent on the modeling of the forward-beamed Cerenkov light and its atmospheric 
scattering. A minimum viewing angle cut is applied to minimize this problem. The cut 
on opening angle between the SDPs is necessary because the geometry obtained by the 
intersection of the SDPs is not well-constrained when the planes are nearly parallel. The 
Gaisser-Hillas profile can be fit to HiRes-1 or HiRes-2 data individually or to both globally. 
While the composition results reported here are from the global fit, the individual fits pro- 
vide additional quality selection criteria. Events for which the per degree of freedom for 
either of the individual fits is larger than 20 are cut. For the global fit, the cut is 15. 

Events for which the disparity between the individual HiRes-1 and HiRes-2 X^a^, fits is 
more than 500 g cm~^ are discarded. Similarly, each HiRes-2 X^^j; from the time binning 
technique described in Section 4.4 is compared to the y^max from the angular binning global 
fit, and events differing by more than 500 g cm~^ are cut. Finally, the geometric uncertainty 
component of the Xmaa; error described in Section 6.3 is required to be smaller than 400 
g cm^^ for the individual HiRes-1 and HiRes-2 fits and smaller than 200 g cm^^ for the 
global fit. Confidence that the fit to Eq. (2) found the correct X^^^ is bolstered when one 
of the sites sees both the rise toward and fall from ^^raax- The bracketing cut ensures this by 
requiring that the measured X^ax is no more than 60 g cm~^ beyond the visible track. 

As described below, in addition to the individual measurements of shower ^^max from 
HiRes-1 and HiRes-2 data, a global y^max fit is performed utilizing all available data. This 
global fit is what is finally used to determine composition. The effect of these cuts on this 
variable is exemplified by Fig. 2 which shows the Monte Carlo global X^oj, resolution as a 
function of ^ and the cut location. Fig. 3 shows the distribution for data after all the 
cuts are applied. 

553 events survive all of the above cuts. The same cuts apphed to the Monte Carlo set 
give a global y^max resolution of 30 g cm~^ and an energy resolution of 13%. (See Figs. 4 and 
5). Fig. 6 shows the distribution of the ratio of log(X^o3;i/X^o3;2) for the data. The width 
of this distribution is a measure of the individual y^ax resolutions. Because of shorter track 
lengths in the single ring, the HiRes-1 y^max resolution is estimated to be 1.4 times worse 
than the HiRes-2 resolution. Taking this into account, and propagating errors, the width of 
the distribution in Fig. 6 corresponds to a HiRes-2 resolution of about 50 g cm~^, which 
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ls in accord with simulations. Fig. 7 shows the same distribution for Monte Carlo events. 
The widths of the distributions are in good agreement. Figs. 8-10 show a representative 
selection of measured event profiles. 

5. EAS SIMULATION 

While cosmic ray hadronic composition presumably can range anywhere between the 
two extremes of pure proton and pure Fe, the 30 g cm~^ resolution of the detector and the 
existence of significant shower fluctuations lead us to compare the data to a simplified two 
component model. Events are generated using CORSIKA 6.005 and 6.010 (Heck 1998), 
using both QGSJetOl (Kalmykov et al. 1997) and SIBYLL 2.1 (Fletcher et al. 1994; Engel 
et al. 1999) hadronic models for both protons and iron nuclei. 

The differences in the hadronic interaction models are evident in the multiplicity, inelas- 
ticity, and hadron-air cross-section they predict. Each of these directly affects the shower de- 
velopment. The QGSJet multiplicity increases as log(E), whereas the multiplicity in SIBYLL 
is below that of QGSJet at relevant energies and rises more slowly than log(E) (Heck et al. 
2001a). 

As a result of this multiplicity dependence, QGSJet showers would be expected to 
develop more quickly than SYBILL events. Both models show an increase of inelasticity 
with energy (Kalmykov et al. 1997; Fletcher et al. 1994; Engel et al. 1999), but QGSJet is 
more inelastic in the UHECR regime (Simpson 2001), again contributing to faster shower 
development (Heck et al. 2001a). However, the hadron-air cross-sections in SIBYLL are 
larger than those of QGSJet at relevant energies. The inelastic p-air cross-section in QGSJet 
rises approximately hnearly with log(E), while in SIBYLL it rises more rapidly (Heck 2001). 
The nucleus-air cross-sections are comparable in magnitude for both models and rise slowly 
with energy (Heck 2001; Simpson 2001). 

In all simulations, the CORSIKA EGS4 option was selected enabhng explicit treat- 
ment of each electromagnetic interaction for particles above a threshold energy. Electrons, 
positrons, and photons were tracked down to energies of 100 keV. Hadrons and muons were 
tracked to 300 MeV. All showers were initiated at 45° to the vertical, with sampling at 5 g 
cm~^ of vertical atmospheric depth, giving bins of about 7 g cm~^ along the shower. 

Because of computational time requirements, simulated UHECR EAS must be generated 
using a "thinning" approximation (Hillas 1997). Numerous studies have shown that setting 
the threshold for thinning at 10~^ of the energy of the primary reduces computation time 
without significantly affecting the results in the mean X^nax and the elongation rate (Simpson 
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2001; Song 2001; Hillas 1997; Pryke 2001). The thinning level for this work was set at 10"^ 

Iron nucleus-initiated showers are expected to have smaller shower-to-shower fluctua- 
tions than proton-initiated showers. Studies with shower generators, including CORSIKA, 
have shown that generating as few as 200 iron showers at a given energy is sufficient to study 
primary composition parameters, whereas 500 proton showers are needed (Pryke 2001). For 
this study, at least 400 iron showers and 500 proton showers were generated using each 
hadronic interaction model in each 0.1 step of log(E/eV) from E — 10^'''^ to 10^° eV. 

The use of thinning enhances fluctuations near shower maximum so that choosing the 
CORSIKA output bin with the largest number of charged particles will often yield a signif- 
icantly incorrect X^ax- The input ^max values used here were obtained using a weighted- 
average smoothing process. 

For this work, a Monte Carlo hbrary containing the complete proflles (number of charged 
particles as a function of atmosphere depth) for all 400-1- iron showers and 500+ proton show- 
ers generated at each energy was constructed. Shower generation in the detector Monte Carlo 
is accomplished by sampling a shower profile in the library, making no a priori assumption 
about the shape of the shower profile. 

6. DETECTOR SIMULATION 

A detailed Monte Carlo simulation of light production, atmospheric transmission, and 
detector response was developed in conjunction with the reconstruction routines, allowing 
the study of the reconstruction code's ability to correctly recover information about the 
primary particle. 

6.1. Monte Ccirlo Program 

The full Monte Carlo simulation of a cosmic ray event begins with the generation of 
an EAS. The longitudinal development of each CORSIKA shower is stored in the library 
described above. After selecting a primary particle energy from the input spectrum, (based 
on the Fly's Eye Stereo Spectrum (Baltrusaitis ct al. 1985)), the Monte Carlo randomly 
selects a shower from the library energy bin closest to the desired energy and interpolates. 
The zenith angle and distance from the detector are then chosen from random distributions. 
In each atmospheric depth bin, scintillation light and Cerenkov light are calculated based on 
the number of charged particles in that bin. The propagation of the hght to the detector is 
then simulated accounting for molecular and aerosol scattering. 
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The photon flux reaching the detector is distributed among the PMTs by a detailed 

ray-tracing program which has been checked by examining the point spread function of stars 
in the night sky (Bergman et al. 2001; Sadowski et al. 2002). The ray tracing accounts 
for PMT cluster obscuration of the mirror, mirror shape, mirror reflectivity, and UV filter 
transmission. Cracks between PMTs are also simulated. A PMT quantum efficiency curve 
(28% at 355 nm), based on specifications provided by the manufacturers, is used to obtain 
the number of photo-electrons. A simulated signal is then generated using the gain of the 
PMT and its pre-amplifier. The triggering conditions are applied, with the resulting output 
stored in the same format as real data. The full timing of each Monte Carlo photon, from 
production at the shower to the PMT face, is stored so that the trigger timing in the output 
is accurate. The Monte Carlo also simulates noise. Sky noise and electronics noise arc added 
to each signal following a Poisson distribution, and random noise tubes arc added to each 
event with the same mean and sigma as in the actual data. The Monte Carlo events are then 
processed by the same stereo and time-binned reconstruction routines used for the data. 



6.2. Atmospheric Comparisons 

Real UHECR events occur in whatever atmospheric conditions present themselves in 
the aperture. Hourly atmospheric parameters are available for most of the data. Fig. 11 
shows the distribution of atmospheric parameters in the database for VAOD of less than or 
equal to 0.1 (corresponding to the data cut). We use this distribution for thrown events. 
The mean scale height is inferred from the mean horizontal attenuation length, the mean 
optical depth, and (Wiencke et al. 2001a) 

AOD = —. (6) 

If no measurement exists in the database, events are reconstructed with the average 
atmospheric description as discussed above. To simulate the error introduced by this, Monte 
Carlo events were generated with input atmospheric parameters sampled from the database 
and then reconstructed with the average parameters. Monte Carlo events were also generated 
with the average atmospheric parameters. Resolutions after quality cuts using the two sets 
were indistinguishable. 
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6.3. Data-Monte Carlo Comparisons 

If the Monte Carlo accurately models the detector, then the application of an event 
selection criterion will have the same effect on Monte Carlo events and data. The Monte 
Carlo can then be used to determine resolution by reconstructing Monte Carlo events and 
comparing the results to the input parameters. Additionally, the effects of selection cuts on 
the resolution can be studied. To determine how well the Monte Carlo models the detector, 
comparisons were made between over 20,000 Monte Carlo events and 926 UHECR events. 

With the cuts described at the beginning of Section 4.5 apphed. Fig. 12 shows the data- 
Monte Carlo comparisons for distribution in Rp. The Monte Carlo histogram in the figure 
represents the results of reconstructing all 20,000+ Monte Carlo events, normahzed to have 
the same area under the curve as the data histogram. The bin-by-bin ratios shown are ratios 
of the data bin height to the normalized Monte Carlo bin height. Fig. 12 is typical of data- 
Monte Carlo comparisons in other variables such as energy, zenith angle, ■?/', tracklength in 
degrees and g cm^^, and maximum single Npg deg~^ m~^ in each event. Good agreement 
between data and Monte Carlo is found in all cases. 

An accurate profile determination depends on an accurate geometry. The SDP-finding 
routine returns the normal to the SDP as well as the uncertainties in each component of 
the normal. For each event, the worst case error in each SDP was propagated through the 
reconstruction to give the uncertainty in Xmax due to the uncertainty in geometry. Fig. 13 
shows the distribution of these uncertainties for X^ax in the data and the Monte Carlo, again 
showing good agreement. 

Since we will be comparing data and Monte Carlo to extract the cosmic ray composi- 
tion, it is vital that the Monte Carlo accurately reproduce the detector and reconstruction 
resolution. This can be demonstrated by examining the differences in the individual detector 
fits. Figs. 14 and 15 show the pull, defined as 2{Xmaxi - ^max2) / (Xmaxi + ^max2), where 
X^axi and Xjnax2 ^re from the individual fits by HiRes-1 and HiRes-2, respectively, after all 
cuts. The nearly-gaussian shape of the pull and the nearly identical pull distribution for 
Monte Carlo and data show that the Monte Carlo resolution represents the real detector 
resolution well. In addition. Figs. 16-18 show the dependence of the data and Monte Carlo 
pull on Xmax, E, Rp and ip angles. There is excellent agreement between data and Monte 
Carlo and no significant dependence of the pull on any of these variables. While the data 
pull distribution cannot be directly used to obtain the global X^ax resolution, the excellent 
agreement between data and MC pulls implies that the MC estimate of the global X^ax 
resolution is reliable. 
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7. ELONGATION RATE RESULT 

The data were binned in energy as shown in Table 2. The statistical errors are the 
standard error of the mean. 

Fig. 19 shows the ER result. The QGSJet and SIBYLL model predictions and the 
HiRcs Prototype result are also indicated. The measured ER is 54.5 ± 6.5 g cm~^ per 
decade (statistical uncertainty only; see Section 7.1), compared to the model predictions of 
50 and 61 g cm^^ per decade for QGSJet protons and iron nuclei, respectively, and 57 and 
59 g cm""^ per decade for SIBYLL protons and iron nuclei, as well as to the HiRes Prototype 
result of 93.0 ± 8.5 (stat) =b 10.5 (sys) g cm~^ per decade. 

7.1. Systematic Uncertainty in Elongation Rate 

Uncertainties in energy do not have a large affect on ER because of the logarithmic 
energy scale. Any systematic uncertainty in X^ax which applies over the entire energy range 
will change the absolute value of X^a^., but will not change the ER. To affect ER, the 
systematic uncertainty must shift X^ai in an energy-dependent way. 

As discussed in Section 4.6, the measured average atmosphere at HiRcs is parameter- 
ized by an aerosol vertical optical depth of 0.04 ± 0.02 (sys). Because the database used to 
obtain this result was used in the reconstruction of about 3/4 of the events, the statistical 
variation about this mean is already represented in the data and the Monte Carlo. The 
effect of the systematic error on the mean is studied by re-processing the data with a dirtier 
atmosphere. Just as in the original processing, the atmospheric database was sampled and 
either the database entry or the standard atmosphere was used, as appropriate. Since hori- 
zontal extinction is measured separately while the scale height is inferred, for each event, the 
aerosol scale height was increased such that the optical depth was increased by 0.02. The 
dominant effect of the dirtier atmosphere was to increase the reconstructed energies. The 
reconstructed X^ax values also decreased slightly, with the two effects combining to steepen 
the ER, as shown in Fig. 20. Because the energy scale is logarithmic the introduced change 
in the ER is small. Since 24% of the data was reconstructed using an average atmosphere, 
we check that deleting this part of the data does not significantly affect our value for the 
elongation rate. This is shown in Fig. 19. 

The modeling of the Cerenkov beam in the reconstruction program can introduce an 
energy-dependent uncertainty in ^max- While the absolute intensity of Cerenkov radiation is 
well known, the effective angular distribution in an EAS depends on the multiple scattering 
of electrons in the atmosphere, which has some uncertainties. To investigate this, the recon- 
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struction code was modified to make the modeled EAS Cerenkov beam 2° wider (consistent 
with one sigma errors in previous measurements), and the data were reprocessed. The ER 
was essentially unchanged. 

8. X^„^ DISTRIBUTION WIDTH RESULT 

Protons are expected to show more shower-to-shower fluctuation than iron nuclei. Fig. 
21 shows that at each of the three energies, both QGSJet and SIBYLL predict that the distri- 
bution of Xmax is wider for proton showers than for iron showers. Thus if the composition is 
changing from Fc to protons as energy increases, for example, then the X^ax distribution at 
lower energies will be significantly narrower than at higher energies. Fig. 22 shows the X„,„a, 
distributions expected from a purely light or a purely heavy fiux with a Fly's Eye Stereo 
spectrum over the entire energy range of interest. For both hadronic interaction models, the 
difference between iron nuclei and protons is clear even over this large energy range. 

Figs. 23-26 show the width of the Xmax distribution. The histograms representing 
the Monte Carlo in Figs. 23-26 were obtained by taking over 4500 showers for each model 
through the complete Monte Carlo and reconstruction routines, subject to the same cuts as 
the data. Nearly 2500 events of each type survived. The areas of the Monte Carlo histograms 
are normalized to the area of the data histogram. 

For Figs. 23 and 24, the data were divided into two energy bins selected such that each 
contained about half of the events. The data shows that the width is not changing with 
energy, indicating that composition is constant or only slowly changing. The width of the 
data distribution in Figs. 25 and 26 indicates that the composition is predominantly light in 
agreement with the elongation rate interpretation. These figures also show the contribution 
of the 24% of the data reconstructed with an average atmosphere. There is no statistically 
significant difference between this sample and the total data sample. 

Figs. 25 and 26 show that the data are consistent with a nearly purely protonic compo- 
sition, especially when compared to the QGSJet model. Assuming a simple two-component 
model where the primary flux is a mixture of protons and iron nuclei, Fig. 27 shows how 
well the model flts the data as a function of fraction of protons. The best fits are at 80% 
protons for QGSJet and 60% for SIBYLL. Fig. 27 also shows the models compared to the 
data. 
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9. SYSTEMATIC UNCERTAINTY IN ABSOLUTE VALUE OF X„„^ 

Fig. 27 suggests that the composition of UHECR is predominantly hght above 10^^ eV. 
However, systematic errors in the absolute value of Xmax could artificially move the measured 
^max values too deep in the atmosphere. The X^ax values for events with energies above 
10^^ eV are of particular interest. Potential contributors to the systematic uncertainty in 
Xmaa; ^re biases introduced by reconstruction, errors in PMT pointing directions, variations 
of the molecular and aerosol component of the atmosphere, and incorrect treatments of the 
Cerenkov beam. 

Resolution studies (Figs. 2 and 3) show that the mean of the reconstructed Xmax 
distribution differs by 5 g cm"^ from the input value. For events with energy greater than 10^^ 
eV, the mean of the distribution shifts by 3 g cm~^. The phototube cluster pointing directions 
have been confirmed by observing stars (Bergman et al. 2001; Sadowski et al. 2002). The 
largest deviation of true pointing direction from the direction used in reconstruction is 0.3°. 
This corresponds to a maximum Xmax error of 15 g cm~^. The molecular component of the 
atmosphere may be different from the U.S. Standard Atmosphere (NOAA 1976) assumed in 
the reconstruction. The variations in the atmosphere as measured by radiosondes launched 
from the Salt Lake City airport has been studied (Martin et al. 1999). In the month of 
the year which varied most from the standard model, the actual pressure differed from 
the model by 8%, which leads to a 10 g cm~^ difference in integrated atmospheric depth. 
Subsequent studies by Y. Fedorova of radiosonde data over several years have confirmed this 
initial study. Using an incorrect atmospheric attenuation could distort the shower profile 
and systematically shift the position of Xmax- Fig. 20 shows that above 10^^ eV, an increase 
as large as 0.02 in VAOD (a 50% change) does not significantly change the mean values of 
^max- As described in Section 7.1, the data were reprocessed with a larger value of Cerenkov 
beam width. The mean difference between the Xmax obtained from reconstruction with the 
standard Cerenkov beam width and the ^max obtained by using a wider beam is less than 1 
g cm~^ and is negligible. 

Table 3 summarizes the systematic uncertainties in X^ax for energies above 10^^ eV. 
Adding the individual uncertainties in quadrature gives an overall systematic uncertainty of 
less than 20 g cm~^. This is less than the detector resolution and much less than the 75 g 
cm~^ difference between proton and iron mean X^ax- 
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10. CONCLUSIONS 

The measured elongation rate result is consistent with a constant or slowly changing 
composition between lO"'^^'^ eV and lO"*^^'^ eV. The data are also in very good agreement with 
the HiRes Prototype data in the region where they nearly overlap. The HiRes Prototype 
result showed a composition change from heavy to hght in the 10^''' to 10^^ eV range, but 
the HiRes data do not show a continuation of this elongation rate. This can be interpreted 
as strong evidence for a transition from a heavy composition to a predominantly light and 
constant or only slowly changing composition above 10^^ eV. The earlier Fly's Eye stereo 
data also shows a transition from a heavy to a light composition though with a broader 
transition region in energy. This is consistent with the present data if systematic errors and 
worse Xjnax resolution quoted in that work are taken into account (Abbasi et al. 2004b). 

In the present study, the widths of the Xmax distributions in the UHECR regime 
strengthen this conclusion. Such a transition is interesting in light of reported structure 
in the UHECR spectrum in this same energy region. Many experiments have seen evidence 
for a second knee in the middle of the 10^^ decade and an ankle structure near 3 x 10^^ eV. 
The change in composition may reflect a change from a dominant galactic CR flux to an 
extragalactic flux which dominates near 10^^ eV. Observation of anisotropy from the galactic 
plane would support this picture. Some evidence for such anisotropy has been reported (Bird 
et al. 1999; Hayashida et al. 1999; BeUido et al. 2001). 

Parts of this analysis and the preparation of this document were performed at Lawrence 
Livermore National Laboratory, which is operated by the University of California for the 
Department of Energy. One of us (Sokolsky) would hke to thank the Guggenheim Foundation 
for support. The HiRes project is supported by the NSF under contract numbers NSF-PHY- 
9322298, NSF-PHY-9904048, NSF-PHY-9974537, NSF-PHY-0098826, NSF-PHY-0305516 
and NSF-PHY-0245428, by the Department of Energy Crant FG03-92ER40732, and by the 
Australian Research Council. The cooperation of Colonel E. Fisher, US Army Dugway 
Proving Ground, and his staff is appreciated. 
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Fig. 1. — The Shower-Detector Plane. The point of the detector and the hne of the shower 
define a plane. 
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Fig. 2. — Global Fit for reconstructed Monte Carlo events. The horizontal line in the top 
plot delineates the cut at — 15- The lower plot shows an expanded view of the region 
remaining after the cut. 
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Fig. 3. — Global Fit for the data. After the application of all cuts, none of the 
a greater than 10. 



data had 
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Fig. 4. — ^max resolution of the detector after all cuts, as determined by the detector 
simulation. X^aj._MC is the thrown X^ax of the shower, and X^aj._recon is the reconstructed 
Y 

^max- 
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Fig. 5. — Energy resolution of the detector after all cuts, as determined by the detector 
simulation. Emc is the thrown energy of the shower, and Erecon is the reconstructed energy. 
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Fig. 6. — Distribution of log(X^aa;i/ Xmax2) in the data. Compare Fig. 7. 
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Fig. 7. — Distribution of log(X^aa.i/ Xmax2) for reconstructed Monte Carlo events. Compare 
Fig. 6. 
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Fig. 8. — Typical reconstructed shower profiles. In the plots on the left, Rpi and Rp2 were 
different by 0.5%. In the plot on the right, they differed by about 5%. The solid curve 
represents the global fit result. 




Fig. 9. — Typical reconstructed shower profiles. In both cases, the track was over 20% 
further from HiRes-1 than HiRes-2. The solid curve represents the global fit result. 
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Fig. 10. — Typical reconstructed shower profiles. Even though HiRes-1 only observed a 
small portion of the shower, the SDP from HiRes-1 stringently constrains the global fit. The 
energy balance depends in detail on atmospheric corrections. In the example on the left, 
a perfect balance was not found, but the measured X^ax was unaffected. The solid curve 
represents the global fit result. 
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Fig. 11. — Distributions of atmospheric parameters. The vertical hnes show the quoted 

average values, which are pulled to the right by measurements which give optical depths 
greater than 0.1 and are therefore not shown. If no aerosols were present, the Horizontal 
Extinction Length would be inifinte (Wiencke et al. 2001a). 
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Fig. 12. — Comparison of reconstructed Rp in the data and the Monte Carlo. The top plots 

show the distributions in the data and the Monte Carlo, and the bottom plots show the 
ratio of data/Monte Carlo in each bin. The similarities of the distributions indicate that 
the thrown events and the Monte Carlo are excellent representations of the data and the 
detector, respectively. 
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Fig. 13. — Magnitude of X^a-r error introduced by uncertainty in the SDP, for both the 
data and MC events, where X„iax-uncert is the X^^^ found after the SDP is shifted by the 
worst-case uncertainty and X^aa;_best is the y^max using the best SDP. The far right bin in 
each plot is an overflow bin. 
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Fig. 14. — Xmax Pull in the data. The pull is defined as 2{Xmaxi - '^max2) / O^maxi + '^max2), 
where X^aa;! and 'y^max2 are from the individual fits by HiRes-1 and HiRes-2, respectively. 
Comparison with Fig. 15 shows excellent agreement between the data and MC pulls. Both 
the pull and log(Xmaa:i / ^max2) are differences between X^axi and y^max2 and are statistically 
equivalent variables, so that Fig. 6 is identical to this plot except for a multiplier on the 
X-axis. 
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Fig. 15. — Xma^ Pull for reconstructed Monte Carlo events. The pull is defined as 2{Xrnaxi 
- ^max2) / O^maxi + ^max2) , wliere Xmaxi and X„iax2 are from the individual fits by HiRes-1 
and HiRes-2, respectively. Comparison with Fig. 14 shows excellent agreement between the 
data and MC pulls. Both the pull and log(X^a3;i / Xmax2) are differences between Xmaxi 
and Xmax2 and are statistically equivalent variables, so that Fig. 7 is identical to this plot 
except for a multiplier on the x-axis. 
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Fig. 16. — Pull as a function of Energy and ^max- The large boxes represent the Monte 
Carlo events, and the small triangles represent the data. 



-39- 



X 



1.5 



0.5 



-0.5 



-1 



a.- ■□; ;q .n''- o ^ g □ □ = = . 



-1.5 



5000 10000 15000 20000 25000 30000 



35000 40000 



1.5 



0.5 



□ .. J.I. r..i...i I i....,f' .1. 1 ■ I, I-. I.I. 1 1 = = ■ ^ ■ ■ 

■ ...t;r..--i.,r-'>-g-. .□■ •□ □. -b □ 

□ • □ n □ n . 



a 



"' •■□.:;%''"t3-(&--n'- • 



>^ -0.5 



□ * □ □. □ ■ 



-1 - 



-1.5 



5000 10000 15000 20000 25000 30000 



35000 40000 
Rp2 (m) 



Fig. 17. — Pull as a function of Rp. The large boxes represent the Monte Carlo events, and 
the small triangles represent the data. 
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Fig. 18. — Pull as a function of ip. The large boxes represent the Monte Carlo events, and 
the small triangles represent the data. 
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Fig. 19. — Elongation rate result. The predictions for QGSJet and SIBYLL protons and 
iron are shown for comparison. The stars show the HiRes Prototype result. The diamonds 
show simulation points calculated by Heck. The best fit to the data and a fit to the 76% of 
the events which have hourly atmospheric corrections are shown. The latter has a slightly 
steeper slope. 
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Fig. 20. — Uncertainty in elongation rate. The circles are obtained by reprocessing the data 
with a dirtier atmoshpere. The x- and y-axes are expanded relative to Fig. 19 to accentuated 
the small difference. The upper two and lower two lines are the model predictions. 
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Fig. 21.— Xmax distributions from CORSIKA for primary UHECR energies of 10^^ ° eV, 
-|^q18.5 g^j^j 10^^° eV. At each energy, the iron distribution is much narrower than the 
proton distribution. 
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Fig. 22. — ^max distributions from CORSIKA for primary UHECR with energies between 
lO^^ o and lO^'^-^ eV. 
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Fig. 23. — Distribution width results, protons. The top plot is for log(E/eV) 18.0 to 18.4 
and the bottom plot is for log(E/eV) 18.4 to 19.4. In each plot, the solid hne is the data, and 
the shaded area represents the 24% of the events reconstructed with the average atmosphere. 
The dashed line is the QGSJet model, and the dotted line is the SIBYLL model. Compare 
Figs. 21 and 24. 
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Fig. 24. — Distribution width results, iron. The top plot is for log(E/eV) 18.0 to 18.4 and 
the bottom plot is for log(E/eV) 18.4 to 19.4. In each plot, the solid hne is the data, and the 

shaded area represents the 24% of the events reconstructed with the average atmosphere. 
The dashed line is the QGSJet model, and the dotted line is the SIBYLL model. Compare 
Figs. 21 and 23. 
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Fig. 25. — All-energy distribution width result, protons. The solid line is the data, and the 
shaded area represents the 24% of the events reconstructed with the average atmosphere. 
The dashed hne is the QGSJet model, and the dotted line is the SIBYLL model. Compare 
Figs. 22 and 26. 
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Fig. 26. — All-energy distribution width result, iron. The solid line is the data, and the 
shaded area represents the 24% of the events reconstructed with the average atmosphere. 
The dashed hne is the QGSJet model, and the dotted line is the SIBYLL model. Compare 
Figs. 22 and 25. 
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Fig. 27. — Simple two-component composition model fit results. On the left, the data points 
show the for a data/model comparison at each per cent proton. On the right, the solid 
line is the data, the shaded area represents the 24% of the events reconstructed with the 
average atmosphere, the dashed line is the QGSJet model with 77% proton, and the dotted 
line is SIBYLL with 57% proton. 
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Table 1. Quality cuts. 



Parameter Cut 



Minimum Viewing Angle (both sites) 10 

Minimum Opening Angle Between SDPs 5 

Individual Site Fit 20 

Global Fit 15 

Timing Fit Xrnax - Global Fit X^„^ 500 

HiRes-1 Fit X^„^ - HiRes-2 Fit 'Krnax 500 

Geometric Uncertainty in Individual Site Fit 400 

Geometric Uncertainty in Global Fit 200 

Bracketing ^rnax ^max within 60 of viewed track 



Note. — Angles are in degrees, are unitless, and X^aa- are in g cm 
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Table 2: Elongation rate data. 



log(E/eV) 


Number of 


Mean 


Mean Xmax 


Bin 


Events 


log(E/eV) 


(g cm 




18.0-18.1 


92 


18.05 ± 


0.01 


711.3 


± 


8.2 


18.1-18.2 


79 


18.15 ± 


0.01 


709.2 


± 


7.5 


18.2-18.3 


82 


18.25 ± 


0.01 


714.5 


± 


9.3 


18.3-18.4 


74 


18.35 ± 


0.01 


729.8 


± 


8.4 


18.4-18.5 


59 


18.45 ± 


0.01 


744.8 


± 


8.5 


18.5-18.6 


34 


18.55 ± 


0.01 


751.1 


± 


11.3 


18.6-18.7 


29 


18.65 ± 


0.01 


734.2 


± 


11.2 


18.7-18.8 


23 


18.74 ± 


0.01 


752.2 


± 


14.1 


18.8-18.9 


23 


18.86 ± 


0.01 


752.3 


± 


8.9 


18.9-19.0 


20 


18.94 ± 


0.01 


757.0 


± 


11.4 


19.0-19.1 


13 


19.05 ± 


0.01 


758.9 


± 


18.0 


19.1-19.2 


14 


19.17 ± 


0.01 


764.4 


± 


7.7 


19.2-19.3 


6 


19.26 ± 


0.01 


781.6 


± 


9.8 


19.3-19.4 


5 


19.35 ± 


0.01 


796.2 


± 


27.7 



Note. — 



Uncertainties are the standard error of the means. 
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Table 3. Systematic uncertainties in X, 



Uncertainty 


g cm ^ 


Pointing Direction 


15 


Atmospheric Variations 


10 


Reconstruction Bias 


5 


Sum in Quadrature 


18.7 



